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ABSTRACT 

This work examines the level of discretization error in simulation-based aerodynamic 
databases and introduces strategies for error control. Simulations are performed using a 
parallel, multi-level Euler solver on embedded-boundary Cartesian meshes. Discretiza 
tion errors in user-selected outputs are estimated using the method of adjoint- weighted 
residuals and we use adaptive mesh refinement to reduce these errors to specified tol- 
erances. Using this framework, we examine the behavior of discretization error 
throughout a token database computed for a NACA 0012 airfoil consisting of 120 
cases. We compare the cost and accuracy of two approaches for aerodynamic database 
generation. In the first approach, mesh adaptation is used to compute all cases in the 
database to a prescribed level of accuracy. The second approach conducts all simula- 
tions using the same computational mesh without adaptation. We quantitatively assess 
the error landscape and computational costs in both databases. This investigation high 
lights sensitivities of the database under a variety of conditions. The presence of tran- 
sonic shocks or the stiffness in the governing equations near the incompressible limit 
are shown to dramatically increase discretization error requiring additional mesh reso- 
lution to control. Results show that such pathologies lead to error levels that vary by 
over factor of 40 when using a fixed mesh throughout the database. Alternatively, con- 
trolling this sensitivity through mesh adaptation leads to mesh sizes which span two 
orders of magnitude. We propose strategies to minimize simulation cost in sensitive 
regions and discuss the role of error-estimation in database quality. 


INTRODUCTION 

The shift to multi-core CPU architectures has rapidly accelerated the growth of su- 
percomputing resources [1], This marked increase in the level of high-performance 
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computing now offers both unprecedented capability and capacity to the aerodynamic 
simulation community [2], These systems are capable of aerodynamic simulations with 
10 8 -10 10 degrees of freedom, offering ever increasing physical fidelity [3,4], While 
such extremely large “capability” simulations are becoming commonplace, the engi- 
neering community has focused on the enormous capacity of these systems through an 
increasing reliance on parametric, trade and optimization studies. In industrial settings, 
the main role of high-end computing is performance database generation. Aerody- 
namic databases with 1 0 3 - 1 0 4 simulations have become common and “production 
CFD” is the mainstay of the workload on this hardware. 

The quality of these databases hinges on the level of the discretization error associ- 
ated with the simulation outputs. Simply put, how much does the mesh influence the 
results? In this work, we use the method of adjoint- weighed residuals to estimate the 
discretization error in each simulation and examine the resulting error distribution. 
Adaptive mesh refinement is used for error control in each case of the database. These 
investigations use the parallel multi-level Cartesian Euler solver developed in refer- 
ences [5] and [6] to produce the aerodynamic data. This simulation package has been 
used extensively on large shared and distributed memory systems and has very good 
parallel scalability [2,3,7], The robustness and automation of this simulation package 
has led to its wide adoption for producing aerodynamic databases in support of engi- 
neering analysis and design [8,9,10,11]. 

The use of adjoint-based error-estimates and adaptive mesh refinement in database 
construction increases the computational cost of each data point because several flow 
and adjoint solutions are required to generate a mesh that satisfies the user-specified 
output tolerances. The adjoint-solver, sensitivity calculation, and error-estimator all use 
the same parallel, multilevel framework as the base Euler solver and all achieve the 
same high level of parallel efficiency [12,13,14], 

In addition to adaptively meshing a simulation, the error-estimation module can 
also be used in a single-pass mode to assess the error in a specific functional on a given 
input mesh. Provided that the Euler simulation on this mesh is sufficiently good, this 
method allows us to accurately assess the discretization-error in a specific output on a 
particular mesh. This approach can be applied to each simulation in an aerodynamic 
database, yielding an a-posteriori estimate of the error landscape for that database. 

Our investigations employ both of these techniques to investigate the role of discre 
tization error in simulation-based aerodynamic data. We begin by briefly reviewing the 
salient features of adjoint-error estimation using a discrete adjoint solver. We then pre- 
sent both fixed-mesh and error-controlled databases for a token geometry in incom- 
pressible, transonic and supersonic flow. Our analysis tracks both error and cost in 
these databases and aims to quantitatively understand the implications of both strate 
gies. 
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Figure 1: Convergence of functional J with mesh refinement showing the definition of total 
error, E (left) and relative error, e (right). 

ERROR ASSESSMENT AND ESTIMATION 

Figure 1 shows two sketches showing convergence of a typical aerodynamic output 
(Cl, Cd, etc.) with mesh refinement. In the sketch at the left of figure l,E is the total 
error in this output due to discretization-error in the numerical solution. This error is 
defined as the difference between the exact functional value,/, and the value obtained 
when evaluating this functional using the discrete solution on the current mesh, /(£/#). 
Rather than attempt to computed directly, we follow the approach of Ref. [16], and 
consider instead the simpler problem of estimating how our discrete evaluation J(Uh) 
would change if we solved on a finer mesh, h. The relative error, e, is sketched at the 
right of figure 1, and is defined as the difference between the functional evaluations on 
the current mesh, H, and the finer, embedded mesh/?. 

e=\J(U H ) -J(Uh)\ (1) 

For a second-order method on a sufficiently smooth solution in the asymptotic 
range, knowing the relative error gives the total error in the output. 

11 4 

E _ e+ -e + — e + --- --e (2) 

Of course, knowledge of the relative error hinges on our ability to evaluate the 
functional on the fine mesh solution, J(Uh). In [12] and [14] we circumvent this diffi- 
culty by approximating the output on the embedded mesh as a function of the coarse 
mesh flow and adjoint solutions. 

J(Uh) * J(US) - (tf ) T R (U«) - W, h - tf) T R (U») (3) 

S v ' V v ' V 

Adjoint Correction Remaining Error 


Where R( ) is the spatial operator of the Euler solver ,xp, is the discrete adjoint so- 
lution and the notation ( )/, is used to indicate prolongation from the coarser to finer 
mesh. 
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References [12] [13] and [15] contain full details of the implementation, and dem- 
onstrate die order of convergence of both the adjoint correction and the remaining er- 
ror terms. In the current work, the remaining error term in eq.(3) is computed by differ- 
encing the linear and quadratic prolongations of the adjoint solution. Earlier publica- 
tions contain detailed verification exercises demonstrating that this fonnulation 
achieves its designed mesh-convergence rate for both the adjoint-correction and the 
remaining-error terms [12, 15]. 


CONSTANT ERROR DATABASE 

As a baseline for investigating the role of discretization error in simulation-based 
aero-data, we start by computing an aerodynamic database to a fixed error-tolerance. 
Our prototype is a simple Mach-a sweep over a NACA 0012 airfoil in 2D. In an effort 
to examine discretization-error and meshing requirements over a wide range of phys- 
ics, our parametric space covers 12 Mach numbers from low subsonic to moderate su- 
personic, Moo = {0. 1-3.0}, and 10 different angles-of-attack in the range {0°-12°}. Fig- 
ure 2 illustrates the extent of this Mach-a space with inset figures (shaded by local 
Mach number) to illustrate flow in various regimes. 

For this study, the objective function was chosen to be a simple unweighted sum of 
lift and drag on the airfoil. 

J(Uh) = Ci+Cj (4) 



Figure 2: Mach-a wind-space covered by a prototypical aero-database showing examples of flow in 
various flight regimes. (12 Mach numbers) x (10 incidence angles) = 120 simulations in the database, 
Mach contours shown for selected cases. 
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Initial Mesh 


980 Cells 



Figure 3: Sample computational meshes in the Mach-a database for NACA 0012 airfoil produced by 
adjoint-based adaptive mesh refinement to control discretization level in functional J= Ci + Cd. 
Initial mesh with 980 cells shown on left. Final meshes had from 10 3 -10 5 cells. 


Aerodynamic Coefficients 


To control discretization-error to a fixed value, every cases in the database was run 
adaptively until the estimate of the remaining-error term in eq.(3) was less than 0.008. 
With the functional in eq.(4), this tolerance translates to 80 counts of drag at zero-lift. 
Figure 3 shows samples of the final meshes required to achieve this tolerance at vari- 
ous points in the database. This error tolerance was purposely chosen relatively loose 
since we hope to be able to achieve it over a wide range of flow conditions. 

Figure 4 shows simulation results for this database through plots of lift, drag and 
moment versus alpha for each Mach number. Reassuringly, these data collapse near the 
incompressible limit, and shock-free inviscid flow shows zero drag for all a. More 
generally, C/ and C m behave linearly where expected, while Q seems to behave nearly 
quadratically. This presentation, however, actually masks the major sensitivities. To see 
this more clearly, figure 5 shows performance landscapes of Ci and Cd now plotted as a 
function of both a and Mach number. While the variation with angle-of-attack is gen- 
erally linear or quadratic, variation with Mach number is much more radical. In this 
view, the steepness of the carpets imply rapid change in the coefficient with either 
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Figure 4: Lift, drag, and quarter-chord pitch-moment for aero-database in figures 2 & 3. 
All cases computed to constant-error tolerance of 0.008 on the approximation of the 
remaining-error in the functional J= Ci + Cd. 
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Figure 5: Performance landscapes of Ci & Cd. Steepness of the carpet indicates output sensitiv- 
ity. All cases computed to constant-error tolerance of 0.008 on the approximation of the 
remaining-error in the functional J=Ci + Cd. 

Mach number or a. The response in either output (C/ or Q) to changes in the flow 
(Mach or alpha) is precisely the definition of output sensitivity. 

Discretization Error and Resolution Requirements 

While figures 4 & 5 show the primary aerodynamic data, figure 6 contains the central 
results for our current investigation of the error landscape. This figure shows quantita- 
tive measures of solution quality and cost. At the left of figure 6, we plot the magnitude 
of the remaining-error term in eq.(3) at each Mach number in the database as a func- 
tion of angle-of-attack. The plot at the right shows the number of cells required to 
achieve this error level also as a function ofM* and a. 

Figure 6 shows the level of remaining error in all the simulations in the database. 
Our error tolerance of 0.008 appears as a shaded (yellow) plane in Mach-a space. In 
examining these data we note that all but one case (M» = 0.5, a = 12°) in the database 
met the error tolerance. Moreover, we see that the adaptive scheme controlled 
discretization-error so that the carpet of remaining error forms a relatively flat surface 
clustered in a narrow band from 0.005 to 0.008. The flatness of this surface is a meas- 
ure of control. 
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Figure 6: (Left) Remaining-error and (right) number-of-cells carpets for constant-error aero- 
database. Shaded plane (yellow) on left shows the error tolerance of 0.008 on functional J = Ci + Cd. 
All but one case achieved the tolerance on remaining-error. 
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Flow solution, M ■*, > 1 Final Mesh 


Adjoint solution, M > 1 



K -R(Uh) = o 
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Figure 7: Confinement of domain in supersonic regime. In supersonic flow, mesh adaptation is con- 
fined to regions behind the bow shock and in front of the limiting characteristic since the remaining 
error term in eq.(3) is zero ahead of the bow shock ( R( (/;,) = 0) and downstream of limiting charac- 
teristic, the output functional has no sensitivity to changes in the flow QPh = 0). 

The carpet plot at the right of figure 6 shows the cell-counts used to achieve this 
level of error. This plot shows a clear division between the subsonic and supersonic 
cases which appears as a sharp step-down in the cell-count carpet located just before 
Mach 1 . The supersonic cases generally achieved the desired level of error with over 
an order of magnitude fewer cells than the subsonic cases. In the subsonic regime, the 
resolution requirements increase with angle-of-attack and we see increasing Mach 
number dependencies approaching the incompressible limit and through the transonic 
regime. 

Figure 7 addresses the disparity in cell-counts between the subsonic and supersonic 
cases. The elliptic nature of subsonic flow means that every point in the domain has the 
potential to influence in our surface functional ./ = Ci + Cd. Since the entire domain 
plays a role, mesh adaptation can extend far upstream or downstream. In supersonic 
flow, however, the domain becomes sharply confined. Upstream of the bow shock, 
Uh = Uoo and since U x satisfies the governing equations, R(f/;,) = 0. Similarly, 
downstream of the limiting characteristic, changes in the flow field cannot effect 
quantities on the body’s surface and thus the adjoint variable is zero. Since the re- 
maining error tenn in eq.(3) is formed by an inner product between the adjoint vari- 
ables and the primal residuals, these zeros confine mesh enrichment to a small re- 
gion of the domain between the bow shock and the limiting characteristic. Similar 
observations were made in ref [14] when discussing figures 9 and 24. This confine- 
ment dramatically limits the problem size in supersonic flow. Moreover, as the Mach 
number increases, this “refinement diamond” contracts in the cross-flow direction and 
the cell count drops with the area change (volume change in three-dimensions). 

In the transonic regime, resolution requirements are driven by the functional’s de- 
pendence on the precise location of the upper surface shock. Since the shock’s position 
is only known to within a cell, accurate shock placement demands fine mesh cells. As 
a result, the simulations at Mach 0.7 are among the most cell- intensive in the database. 
By Mach 0.9 the shocks are attached to the trailing-edge, shock location is unambigu- 
ous, and these cases are easy by comparison. 
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Figure 8: Final meshes and isobars in discrete solutions required to achieve an error- 
tolerance of 0.008 on functionalJ= Cl + Coin near-incompressible flow. A // = {0.1, 0.2, 

0.3} and a = 4°. 

While cell-counts in the trans- and supersonic regimes are driven by flow physics, 
near the incompressible limit they are driven by our model of the governing equations. 
These simulations were preformed using the compressible form of the governing equa- 
tions without any low-Mach preconditioning. Near the incompressible limit, the pres- 
sure coefficient becomes independent of Mach number giving rise to the well-known 
1/M 2 incompressible scaling. Figure 8 shows isobars in the discrete solution for the 
three lowest Mach numbers. Isobar levels at M» of 0.1, 0.2 and 0.3 have been chosen 
to illustrate the approach to self-similarity in these nearly incompressible flows. Pres- 
sure signals at Moo = 0.1 are 9 times weaker than their counterparts at M» = 0.3 in ac- 
cordance with the inverse Mach-squared scaling. This represents nearly a 10 fold de- 
crease in signal strength and therefore meeting the same error-tolerance on the aerody- 
namic coefficients requires correspondingly higher mesh resolution. This stiffness is 
simply an artifact of solving the un-preconditioned governing equations and the fact 
that our functional is based upon coefficients which depend onq^. 
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Figure 9: Computational mesh with -7000 
cells used for all cases in the fixed-mesh data- 
base. 



Figure 10: Drag in database computed on fix 
mesh. Inset shows poor prediction of drag at 
Mach near a = 0. 


CONSTANT MESH DATABASE 

In contrast to the adaptively refined meshes in the preceding section, we now ex- 
amine a database computed on a fixed computational mesh. Figure 9 shows a mesh 
with -7000 control volumes. This mesh is similar in cell-density and cell-count to 
some of the adaptive meshes found at moderate subsonic Mach numbers (see fig 3.) 
and is representative of a “best practices” mesh for this airfoil. Using the same mesh 
for all simulations in a database is common practice in production CFD making it im- 
portant to understand the implications of this approach in terms of data quality and dis- 
cretization error. This fixed-mesh database used the same computational mesh for all 
120 simulations in the database. Figure 10 shows the prediction of drag across the en- 
tire database as a function of angle-of-attack for all Mach numbers. In comparing these 
results with those in Fig. 4 we notice immediately the poor prediction in the subsonic 
regime at low incidence angles. Spurious drag in this region increases as the flow be- 
comes more incompressible, and increases with a. Given the earlier discussion of flow 
sensitivity in this regime, these results are hardly surprising, but the quantification is 
illustrative. In shock-free flow, inviscid flow theory predicts zero drag for this airfoil. 
Instead, Fig. 10 shows finite drag for all these cases - with results nicely sorted with the 
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Figure 11: Performance landscapes of C; & O/ for database computed on a fixed mesh with -7000 
control volumes. Note poor prediction of drag in subsonic shock-free flow. Symbol color indicates 
level of remaining-error in simulations relative to tolerance of 0.008 on on functionalJ= Cl + Co. 
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inverse of Mach num- 
ber, as we approach 
the incompressible 
limit. Figure 11 shows 
performance land- 
scapes with Ci and Cd 
plotted as functions of 
Mach and angle-of- 
attack. Its interesting 
to contrast these with 
their counterparts in 
fig. 5 and the poor 
prediction of drag in 
subsonic shock-free 
flow is obvious. Sym- 
bols on the plot in fig- 
ure 5 are colored by 
the value of 
remaining-error in the 
simulations. Red symbols mark cases with error at least 50% in excess of our 0.008 
tolerance on J= Ci + Cd, and virtually all of the subsonic cases are red. 

Figure 12 gives a more quantitative understanding of the error-landscape for the 
constant-mesh database. To facilitate comparison, these data are plotted using the same 
scale as the error-carpet in the error-controlled database (fig. 6). The yellow plane in 
fig. 12 shows the tolerance of 0.008 on J. The undersolved simulations (red symbols) 
are above this plane, while the oversolved ones (blue symbols) lie below. White sym- 
bols designate cases that lie within 50% of our desired error tolerance. In general the 
cases meet or exceed the error tolerance at Mach 0.9 and above. Although all cases 
were computed on the same mesh, the remaining-error data in fig. 12 vary from 0.016 
to 0.06. This means that there is a factor of 40 difference in accuracy of our lift and 
drag functional between the best and worst cases in the database. In general we ob- 
serve that the subsonic cases have 5-10 times more error than the supersonic cases, and 
that the errors are worse in regions of the domain where we the governing equations 
are stiff (transonic regime and near the incompressible limit). 

This last observation underscores the duality between meshing requirements and 
discretization error. Figure 13 clarifies further. On the left is the “meshing require- 
ments” carpet plot from the error-controlled database replotted from figure 5. On the 
right, we replot the error carpet from the fixed-mesh database in figure 12, only this 
time, using a logarithmic scale on the vertical axis. Comparing the two reveals remark- 
able similarities. We see that the cases with the highest error on the fixed mesh are ex- 
actly the same as those which the adjoint has identified as having the most stringent 
meshing requirements. Cases with high sensitivity have both high meshing require- 
ments and high-error when the fixed-mesh is used for their simulations. Of course, this 
is not surprising, but we can now quantify the error range and realize that even for a 



Mach Number 

Figure 12: Remaining error landscape of fixed-mesh database. Yellow 
plane indicates error tolerance of 0.008 in functional J= Cl + Cd. 
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Meshing Requirements Error Landscape 

(Constant-Error Database) (Constant-mesh Database) 



Figure 13: Duality between meshing requirements in constant-mesh database and remaining error 
landscape of constant-mesh database (logio scale). 

simple 2D airfoil, using a fixed mesh results in output error that spans nearly 2 orders 
of magnitude. Conversely, attempting to control this error, even in 2D, requires meshes 
with with cell-counts that vary over two orders of magnitude (1600-200000 cells). 

As discussed earlier, cell-counts near the incompressible limit are driven by the 
inverse Mach-squared incompressible scaling. Cell-counts in the transonic range are 
driven by the need to resolve the precise location of the shocks, and when this position 
is only known to within one cell, //-refinement is the only meaningful way of improv- 
ing it. In the high-Mach cases, the domain is substantially confined dramatically reduc- 
ing meshing requirements, but a certain number of cells are still required to resolve the 
shock and flowfield behind the shock. 


ALTERNATIVE STRATEGIES 

In many ways, blindly insisting that all cases in a database achieve a constant error 
is hardly better than blindly applying the same mesh to all cases. Given the range of 
cell-counts required to achieve constant error, the hardest few percent of cases will 
completely dominate the computational cost of the database. In all likelihood a data- 
base which truly has uniform error may not be demanded for engineering use. Trajec- 
tory simulations, for example have their own set of sensitivities. Cases where the vehi- 
cle is likely to spend the majority of its time generally require tighter error tolerances 
than comers of the flight envelope that may never be reached. Under these circum- 
stances, simply capping the maximum number of cells may be a useful technique for 
avoiding the needless expense chasing precision in comers of the flight envelope 
which are unimportant in trajectory simulations or for Guidance and Control system 
development. On the other hand, its also interesting to consider a relative tolerance 
where the functional’s value is linked to its error tolerance. 
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CONCLUSIONS AND FUTURE WORK 

In order to examine the level of discretization error in simulation-based aerody- 
namic data, this work compared error levels in two aerodynamic databases using a par- 
allel, multi-level Euler solver on embedded-boundary Cartesian meshes. Discretization 
errors in user-selected outputs were estimated using the method of adjoint-weighted 
residuals and adaptive mesh refinement was employed to reduce these errors to speci- 
fied tolerances. Using this framework, the investigations examined the behavior of dis- 
cretization error throughout a token database computed for a NACA 0012 airfoil con- 
sisting of 120 cases. The accuracy and computational expense were compared for two 
approaches for aerodynamic database generation. The first approach used mesh adapta 
tion to produce an “error controlled” database where all cases in the database achieved 
a prescribed level of accuracy. The second approach conducted all simulations using 
the same computational mesh without adaptation. Error landscapes in both approaches 
were quantitatively assessed along with the computational costs of both approaches. 

This investigation highlighted the interplay between flow field sensitivity, discreti- 
zation error and meshing requirements. Aerodynamic coefficients in the error- 
controlled databases demonstrated textbook behavior. However, it was also shown that 
due to varying flow sensitivity, maintaining a constant error tolerance can require 
widely disparate mesh sizes. The presence of transonic shocks or the stiffness in the 
governing equations near the incompressible limit were shown to dramatically increase 
discretization error requiring additional mesh resolution to control. Results show that 
such pathologies lead to error levels that vary by over factor of 40 when using a fixed 
mesh throughout the database. Alternatively, controlling this sensitivity through mesh 
adaptation leads to mesh sizes which span two orders of magnitude. 

We intend to repeat this study in 3D to understand and quantify the impact of di- 
mensionality on both error levels and meshing requirements. Furthermore we intend to 
examine the usefulness of capping cell-counts to control database cost and examine the 
utility of relative error-tolerances as a middle-ground between accuracy and expense. 
Clearly quantitative error estimates on simulation data are of great value and error- 
estimates are reasonable expectation using techniques that are already available. The 
dual nature of meshing requirements and discretization error implies that even in the 
absence of mesh adaptation, the error estimates offered by the method of adjoint- 
weighted residuals adjoint provides insight into both flow sensitivity and meshing re- 
quirements. 
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